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ABSTRACT 

Using the hydrodynamic code ZEUS, we perform 2D simulations to determine 
the fate of the gas ejected by massive stars within super star clusters. It turns out 
that the outcome depends mainly on the mass and radius of the cluster. In the 
case of less massive clusters, a hot high velocity (~ 1000 km s -1 ) stationary wind 
develops and the metals injected by supernovae are dispersed to large distances 
from the cluster. On the other hand, the density of the thermalized ejecta within 
massive and compact clusters is sufficiently large as to immediately provoke the 
onset of thermal instabilities. These deplete, particularly in the central densest 
regions, the pressure and the pressure gradient required to establish a stationary 
wind, and instead the thermally unstable parcels of gas are rapidly compressed, 
by a plethora of re-pressurizing shocks, into compact high density condensations. 
Most of these are unable to leave the cluster volume and thus accumulate to 
eventually feed further generations of star formation. 
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The simulations cover an important fraction of the parameter-space, which 
allows us to estimate the fraction of the reinserted gas which accumulates within 
the cluster and the fraction that leaves the cluster as a function of the cluster 
mechanical luminosity, the cluster size and heating efficiency. 

Subject headings: Galaxies: star clusters — ISM: bubbles — ISM: HII regions - 
ISM 



1. Introduction 



Due to their stellar mass, which in some cases amounts to several million M , and their 
compactness, as they only span a few parsecs, young (< 10 Myr) super star clusters (SSCs) 
represent the most spectacular and dominant mode of star formation in starburst and in 



ter ac ting galaxies (jO'Connell et al.lll995l : lHolll997l : IWhitmore &: Schweizerlll995l : iMelo et al. 



20051 ). SSCs have been detected in the optical, UV and X-rays, and also in the radio contin- 
uum and IR regimes, as some of them are d eeply embedded behind dense obscuring material, 
leading to powerful ultra-dense HII regions (jKobulnicky fc Johnsonl|l999l : iGilbert &: Graham 
2007h . 



On theoretical grounds, it has been inferred that such extreme modes of star formation 
should lead to several tens of thousands of massive stars, all of them known to rapidly 
(< 50 Myr) reinsert, through stellar winds and supernova (SN) explosions, a large fraction 
of their mass back into the ISM. The first (a diabatic) approach t o the hydrodynamics of 
the matter reinserted within a SSC is due to I Chevalier fc Clegs (119851 . hereafter CC85). 
In their model, the stellar sources of mass and energy, are assumed to be equally spaced 
within the SSC volume of radius R$c- They also assumed that all of the kinetic energy 
provided by massive stars is immediately, and in situ, thermalized via random collisions 
of the ejecta from neighboring sources. This results in energy and mass deposition rate 
densities q e = (3L sc )/(47ri?g C ) and q m = (3Msc)/(4vri?g C ), respectively, where L sc and Msc 
are the cluster mechanical luminosity and mass deposition rates. These assumptions lead 
to a high temperatures gas (T > 10 7 K) in which the interstellar cooling law is close to 
its minimum value, and this justifies their adiabatic assumption. In the adiabatic model of 
CC85, the thermalized hot gas rapidly settles into almost constant density and temperature 
distributions, although a slight outward pressure gradient establishes a particular velocity 
distribution with the stagnation point (i.e. zero velocity) at the cluster center. The gas 
velocity increases then almost linearly with radius to reach the sound speed (esc) right at 
the cluster surface and then streams into the surrounding lower pressure ambient medium to 
reach its terminal velocity (vaoo ~ 2csc) ; while the density and temperature of the outflow, 
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the cluster wind, decrease as r~ 2 and r~ 4//3 , respectively. The solution of such a stationary 
outflow depends on three variables: the cluster size (Rsc), the mass deposition rate (Msc) 
and the mechanical luminosity of the cluster (Lsc)- Knowledge of these three variables 
allows one to solve the hydrodynamic equations analytically and workout the run of density, 
temperature and velocity of the stationary outflow. Note that Msc is usually replaced by 
the adiabatic terminal speed (uaoo — (2Lsc/Msc 



The model then yields a stationary flow in which the matter reinserted by the evolving 
massive stars (Msc) equals the amount of matter feeding the cluster wind (47ri?g C psc c sc); 
where psc is the reinserted gas density value at the star cluster surface. As Lsc and Msc 
increase linearly with the cluster mass (Lsc ~ Msc, Msc ~ Msc), the adiabatic model 
predicts that the more massive a cluster is, the more powerful is its res ultant wind. Thi s 
latter conclusion breaks down if one relaxes the adiabatic assumption (see lSilich et al.ll2004j ). 
More massive clusters deposit larger amounts of matter and thereby deliver a sufficiently 
high density, psc = ^sc/(4vri?g C ), to provoke strong radiative cooling. Thus, as radiative 
cooling (Q = n 2 A(T, Z); where A(T, Z) is the cooling function) is proportional to M| c and 
Lsc is proportional to Msc, there is a threshold mechanical luminosity (for given i?sc and 
M S c), above which strong radiative cooling takes over despite the large temperatures and 
the minimum value of the interstellar cooling law at these temperatures. 

Note that details of the thermaliz ation process have b e en largely ignored although more 
recent formulations of the problem (jWunsch et al.l 120071 ; ISilich et al.l 120071 ) have inferred 
that a significant fraction of the deposited mechanical energy could be radiated away as 
soon as it is inserted. In such cases, only a fraction of the deposited energy remains avail- 
able to heat the thermalized matter. This fraction is called the h eating efficiency, 77 . It is 
assumed by different authors to hav e values between 0.01 and 1 (IBradamante et al.l Il998 



sive compact SSCs (see ISilich et al. 



Melioli fc de Gouveia Dal Pinoll2004j). an d shown to acquire small values in the case of mas- 



20071 ). 



Figured] presents the threshold mechanical luminosity found by lSilich et al.l (120041 ). here 
re-calculated for three different values of the heating efficiency 77. Clusters with a mechanical 
luminosity (or mass) far below this line evolve in the quasi-adiabatic regime. For these, the 
Chevalier & Clegg mode l provides a good approximation to the structure and hydrodynamics 
of the star cluster wind (jCanto et al.ll2000l ; iRaga et al.ll200ll ). Figure [2] displays radial profiles 
of temperature, particle density, pressure and velocity, typical of such steady state winds. 

For clusters with a mechanical luminosity close to the threshold value, the temperature 
distribution with in the stationary wind becomes very different from that predicted by the 
adiabatic model (ISilich et al.ll2004l ). This is because, as soon as the temperature of the wind 
decreases to approximately ~ 10 6 K, radiative cooling begins to increase sharply (mainly due 
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to f-b and b-b transitions of heavy elements) and the temperature in the free wind region 
drops rapidly several orders of magnitude. Nevertheless such clusters manage to eject all the 
matter deposited by stellar winds and supernovae, and thereby sustain a stationary wind. 

For clusters above the threshold line, radiative cooling becomes a dominant factor. 
Radiative cooling affects first the central densest regions causing a sudden drop in pressure. 
This promotes the shift of th e stagnation point out of the cluster center. Such clusters 
adhere to a bimodal solution (ITenorio-Tagle et al.l 120071 ; IWiinsch et al.l 120071 ) in which the 
stagnation radius (R s t) defines two different regions within the cluster volume (see Figure [3]). 
On one hand, there is an outer shell in which the deposited energy, although affected by 
strong radiative cooling, is still able to drive an outward stationary wind, by making the gas 
velocity acquire the sound speed exactly at the cluster surface. On the other hand, the matter 
deposited inside the stagnation radius is strongly affected by radiative cooling and becomes 
thermally unstable. The instability leads to a rapid and continuous loss of energy from large 
and small parcels of gas, thereby reducing their temperature and pressure. These events 
lead immediately to the formation of strong shocks that emanate from the hot, high pressure 
regions and are driven into the cold parcels of gas in order to restore their pressure. These 
have been termed re-pressurizing shocks, and have been i nvoked in several as trophysical 
circumstances, such as the formation of globular clusters (IVietri &: Pescd Il995l ). and also 
in the cooling of super nova matter within s uperbubbles, leading to highly metallic droplets 
falling onto the galaxy (jTenorio-Tagle!ll996l ) . In the context of the matter cooling inside the 
superstar cluster stagnation radius, the re -pressurizing shocks have b een first modelled by 
means of ID numerical hydrodynamics in ITenorio-Tagle et al.l (120071 ). The re-pressurizing 
shocks rapidly compress the cold gas, enhancing its density while reducing its volume, until 
the cold condensations acquire again the thermal pressure value of the hot gas. Given the 
initially similar values of density ahead of and behind the shocks, one can show that their 
velocity is only a function of the temperature T of the hot gas (Vrp ~ [(kT) / (fj,m p )} ' 5 , 
where k is the Boltzmann constant, fim p is the mean gas-particle mass, and m p is the proton 
mass), and thus cold parcels of gas within the SSC are rapidly driven into small high density 
condensations. 



The continuous occurrence of thermal instabilities results in the accumulation of mass 
in this region and in a very chaotic, highly non-stationary hydrodynamical pattern, with a 
number of radiative shocks and cooling fronts propagating within the cluster volume. The 
continuous accumulation of ther mally unstable matter mus t finally lead to its re-processing 
into further generations of stars (ITenorio-Tagle et al.l 120051 ). Unfortunately, ID simulations 
do not allow one to reach an adequate understanding of the physics within massive clusters 
in the bimodal regime, nor to make realistic predictions regarding the fate of the matter 
reinserted by massive stars. In order to develop a more realistic model, here we present 
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detailed 2D hydrodynamic simulations of the gaseous flows that result from the thermaliza- 
tion of the kinetic energy supplied by stellar winds and supernovae ejecta inside the volume 
occupied by the stellar cluster. We focus on clusters evolving in the bimodal regime. The 
major result from these simulations is that the central zones of clusters evolving in the bi- 
modal regime accumulate large amounts of matter in the form of warm (T ~ 10 4 K) high 
density condensations embedded into a hot plasma of much lower density. The amount 
of accumulated matter depends on the excess star cluster mechanical luminosity over the 
threshold value and grows with time, unless the accumulation of reinserted matter inside the 
stagnation radius is compensated by secondary star formation. Our results also show that 
the stagnation surface itself has a very complicated dynamic morphology that continuously 
changes with time. Nevertheless, the average radius of the stagnation surface remains close 
to that predicted by ID and semi-analytic calculations. 

The paper is organized as follows. Section 2 describes the numerical model and discusses 
the input physics. Section 3 present the results from our numerical simulations and compares 
them with the semi-analytical results and previous ID simulations. In section 4 we discuss 
our results and section 5 gives a summary of our findings. 
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Fig. 1. — The threshold lines. Clusters above the threshold line evolve in a bimodal regime 
in which the injected matter is accumulated inside their densest inner regions while the outer 
zones develop a strongly radiative stationary wind. The threshold lines were calculated for 
the adiabatic terminal velocity Vaoo = 1000 km s _1 and three different heating efficiencies 
r) = 1.0,0.3 and 0.1 denoted with solid, dashed and dotted lines, respectively. The 2D 
simulations (see Table [1]) are represented by symbols. Different symbols denote different 
heating efficiencies (and thus are to be compared with the corresponding threshold line) 
7] = 1.0 (triangles), r] = 0.3 (circles), and r\ = 0.1 (plus signs). The secondary y-axis shows 
the approximate mass of the cluster obtained using a relation M sc = (L sc /3 x f0 40 )10 6 M 
rtLeitherer et alJll999h . 
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Fig. 2. — Model 1 (L sc /L crit = 0.5): radial profiles of the wind particle density, n w (note the 
units are m~ 3 to fit within the figure), temperature (T w ), pressure (P,„) and radial velocity 
(v w ). The thin lines represent the semi-analytical solution (ISilich et al.ll2004j ). the thick lines 
are results from the 2D simulation (model 1) at t = 0.25 Myr. The simulation, after a short 
initial relaxation period, stays perfectly stationary and spherically symmetric at all times. 
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Fig. 3. — The internal structure of SSCs in the bimodal regime. The cluster volume presents 
two distinct regions. Below the stagnation radius R st : is the thermally unstable region 
and above it lies the outer quasi-stationary wind region. The radial profiles of the wind 
particle density, temperature, pressure and radial velocity were obtained as the axial time 
averages values of the 2D simulation (model 3) between 0.4 and 0.8 Myr. Only the hot gas 
(T > 2 x 10 4 ) was taken into consideration for this plot. The values of temperature and 
pressure at the stagnation radius, obtained from the semi-analytical model (thin solid lines 
across the inner region), are given with horizontal lines. 
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2. The numerical approach 



The numerical models pres ented here are based on the finite difference Eulerian hydro- 
dynamic code ZEUD3D v3.4.2 Jstone & Normanlll992h . All simulations have been carried 
out on a spherical (r, 9) grid, symmetric along the 0-coordinate. We have set the radial size 
of grid cells Ar proportional to the radial coordinate r which ensures that all grid-cells have 
Ar ~ rA6. Another advantage of this radially scaled grid is that the resolution is higher at 
smaller radii (inside the cluster) where thermal instabilities are expected. 

The cooling routine accounts for extremely fast cooling both in the wind or within the 
SSC volume. The change of internal energy e due to cooling is 

de 



dt 



-niTi e A(T, Z) 



1 



cool 



where rti and n e are ion and electron densities, respectively. We compute them from the 



n, 



pj ■ 



where ii\ ^m v is the average ion mass. We assume 



gas density p as rii 

[i = 0.609 and /Xi on = 1.27 for all computed models. A(T, Z) is the cooling function. We 
use Raymon d fc Co x cooli ng function which has been supplemented with new elements and 
tabulated bv lPlewal Jl995h . 



The RHS of equation (CQ) is evaluated in the middle of time-steps to maintain the second 
order accuracy of the code and the energy conservation equation is then solved iteratively 
using the Brendt algorithm which is more stable and accurate than the Newton-Raphson 
method originally used in ZEUS. 

The cooling rate has to be included in the computation of the time-step, otherwise rapid 
cooling not resolved in time may lead to the occurrence of negative temperatures. A common 
way of solving this is to limit the amount of internal energy which can be radiated away 
during one time-step by setting 

e 



dt, 



cool 



de I 

dt I cool 



where e is a safety factor smaller than unity (see e.g. ISuttner et al.l 119971 ; where e 
used). In this work we use e = 0.25. 



(2) 
0.3 was 



dt 



The global time-step, dt, is computed as follows: 
dt = dtftD, for <it coo i > dtim; 

dt = dt cooh for dt UD > dt cooi >5x dt HD ; (3) 
dt = 5 x g^hd, for 5 x dt^D > dt coo \ (local substeps dt sn \ } = dt coo \); 

where g?£hd is the " hydrodynamic" time-step resulting from the Courant-Friedrich-Levi cri- 
terion and S is the minimum fraction of g^hd to which the global time-step dt is allowed to 
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drop. If, due to the cooling condition, a certain cell requires an even smaller time-step (i.e. 
dt cool < 5^hd)i the energy equation is integrated numerically in that cell, using rft su b = dt coo i 
but assuming that during this time the density and temperature in the cell are not substan- 
tially affected by interactions with neighboring cells. This ensures that CPU time is not 
wasted in cells where a high time resolution is not required. To determine a reasonable value 
of 5 we ran several tests on a low resolution grid (150 x 56) and found that there are no 
appreciable differences for 5 < 0.3. Therefore, we use 5 = 0.1. 

In order to simulate the effect of the stellar UV radiation field, in some of the simulations 
(see Table [I]), we do not allow the gas temperature to drop below T\ [m = 10 4 K. This is 
equivalent to assuming that there are sufficient UV photons to ionize the dense thermally 
unstable matter which otherwise would cool to much lower temperatures. 

The wind source was modelled by a continuous replenishment of mass and internal 
energy in all cells within the cluster volume, at rates q m = (3Msc)/(47ri?g C ) and q e = 
(3Lsc)/(47r_Rg C ), respectively. The procedure applied to each cell within the cluster volume 
at every time-step is: 

1. the density and the total energy in a given cell are saved to p Q id and e to t,oid 

2. the new mass is inserted p new = p id + (1 + A noise ()q m dt, the velocity is corrected so 
that the momentum is conserved v new = v i d p id/ Pnew 

3. the internal energy is corrected to conserve the total energy ei im id = e to t,oid _ PnewVn ew /2 

4. the new energy is inserted in a form of internal energy e i)new = e i)mi d + (1 + A noise ()q e dt 

where ( is a pseudo-random number from the interval ( — 1,1) generated each time it is used, 
and A n oi S e is the relative amplitude of the noise. The inclusion of a noise term is necessary 
to break the artificial spherical symmetry imposed by the initial conditions (see below). The 
model is very robust with respect to A noise . Test runs with v4 no i se = 0.01, 0.05, 0.1 and 0.5 lead 
to very similar general properties (mass fluxes at boundaries, number of fragments formed 
and their approximate sizes) in all models. The only noticeable difference is the duration of 
the initial relaxation period required to break the initial symmetry. We use A no i se = 0.1 for 
all models described in this paper. 

The boundary conditions are set open at both r-boundaries and periodic at both 9- 
boundaries. The open inner r-boundary allows the dense clumps that are not ejected from the 
cluster (through its boundary at Rsc) to leave the computational domain. Otherwise, they 
would accumulate within the cluster and eventually fill its whole volume; this happens for 
some models with high Lsc and low rj, even with the open inner r-boundary (see Section [375]) . 
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lt is unphysical, because the accumulated mass exceeds the amount which can be ionized 
by the available UV photons; therefore it should cool down, become gravitationally unstable 
and collapse into stars. Moreover, the accumulated mass ultimately becomes so high that 
it would be gravitationally unstable even at T = 10 4 K. However, determining the fate of 
the dense mass properly would need a model of star formation which takes into account 
radiation transfer and self-gravity: this is not included in our current numerical model. 

Two different initial conditions are used for the two series of models presented here (see 
Tabled]). Models 1-6 and 12 with R$c = 10 pc h ave as initial cond ition the spherically 



symmetric semi-analytical solution with Lsc = -^crit (ISilich et al.ll2004j ). Models 7-11 and 
13, with Rsc = 3 pc, use as initial condition the semi-analytical solution with the appropriate 
L$c- As this is only fully defined for r > i? st , at t — the central region r < R st is filled 
with a zero velocity, constant density and temperature gas, with values equal to those at the 
stagnation radius p = p(R s t) and T = T(i? st ), respectively. The advantage of this approach 
is that such conditions are closer to the quasi-stationary state and therefore it takes a shorter 
time to reach it. 



3. Results 

We have calculated two series of models (see Tabled]). Models in the first series all have a 
cluster radius Rsc — 10 pc and an assumed heating efficiency rj = 1. The cluster mechanical 
luminosities considered are: 2.5 x 10 41 , 10 42 , 10 43 and 10 44 erg s _1 , which result in values of 
L sc /L crit = 0.5, 2, 20 and 200, respectively. The model with L sc /L crit = 2 was computed 
with three different numerical resolutions 150 x 56, 300 x 112 and 600 x 224 to check and 
establish convergence. The computational domain extents radially from R^ = 2 pc (the 
inner boundary) to i?oB = 30 pc (the outer boundary) and from 6*lb = vr/2 — 0.5 (left 
boundary) to ^rb = 7r/2 + 0.5 (right boundary) in the axial direction. 

For the second series of models, we assume a more compact cluster and more realistic 
cluster parameters by setting R sc = 3 pc and r\ = 0.1 or 0.3. We ran 6 models with 3 
different ratios Lsc/^crit = 2-5, 25 and 250. In these cases the radial extent of the grid goes 
from Rib = 0.5 pc to Rob = 10 pc, and the axial extent goes from n/3 to 27r/3. 

In all cases we assumed an adiabatic terminal velocity is t>A,oo = yj~fif^ = 1000 km/s, 
and the lower temperature limit was set equal to T\ im = 10 4 K, or 10 2 K (models 12 and 
13). In cases with an r] < 1, t>A,oo, was replaced by t>A,oo = \J 2 }f^ ■ Note that in all cases, 
radiative cooling lowers the outflow velocity at the grid outer boundary to somewhat smaller 
values. 
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No. 


Lsc 
[erg/s] 


Rsc 

[pc] 


V 


Lsc 

^crit 


grid 


T 

± lim 


[PC] 


T st 
[10 6 K] 


M ou t 


T ■ 

± mm 


1 


2.5 x 10 41 


10 


1 


0.5 


300 x 112 


10 4 





10.4 


1.00 


10 4 


2 


10 42 


10 


1 


2 


150 x 56 


10 4 


6.056 


10.4 


0.81 


10 4 


3 


10 42 


10 


1 


2 


300 x 112 


10 4 


6.056 


10.4 


0.84 


10 4 


4 


10 42 


10 


1 


2 


600 x 224 


10 4 


6.056 


10.4 


0.79 


10 4 


5 


10 43 


10 


1 


20 


300 x 112 


10 4 


8.991 


10.4 


0.39 


10 4 


6 


10 44 


10 


1 


200 


300 x 112 


10 4 


9.696 


10.4 


0.19 


10 4 


7 


10 39 


3 


0.1 


2.5 


300 x 104 


10 4 


1.872 


1.10 


0.77 


10 4 


8 


10 40 


3 


0.3 


2.5 


300 x 104 


10 4 


1.860 


3.48 


0.84 


10 4 


9 


10 40 


3 


0.1 


25 


300 x 104 


10 4 


2.709 


1.10 


0.98 


10 4 


10 


10 41 


3 


0.3 


25 


300 x 104 


10 4 


2.707 


3.48 


0.83 


10 4 


11 


10 41 


3 


0.1 


250 


300 x 104 


10 4 


2.912 


1.10 


1.00 


10 4 


12 


10 43 


10 


1 


20 


300 x 112 


10 2 


8.991 


10.4 


0.38 


10 2 


13 


10 40 


3 


0.1 


25 


300 x 104 


10 2 


2.709 


1.10 


0.44 


10 2 



Table 1: The set of computed models. The values of R st and T st were determined by means 
of semi-analytic models, the values of M out /Msc are the mass flux through the cluster border 
in the 2D simulations, averaged over the time interval 0.1 — 0.8 Myr for Rsc = 10 pc models 
and 1 — 2 Myr for i? S c = 3 pc models. 
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3.1. Model 1, L sc /L crit = 0.5 

This model was computed in order to test the numerical code against the semi-analytical 
solution which is known for Lsc < L cr it. Despite the perturbations of the deposited mass 
and energy, the flow is perfectly stationary. The resultant radial density, temperature and 
velocity profiles are shown in Figure [2] where they are compared to the semi- analytic solution. 
The agreement is very good, despite the fact that the 2-D model does not calculate the central 
sphere of radius Rib, and this induces a small discrepancy in the velocity in the inner cluster 
regions. 

3.2. Model 5, L sc /L crit = 20 

We continue with a detailed description of model 5 which exhibits the typical hydro- 
dynamic behavior for clusters above the threshold line, Lsc > L cr - lt . The model starts with 
the semi-analytical solution for Lsc — L cr i t and from then onwards, at each time step, the 
mass and energy input rates are consistent with the selected ratio L sc /L crit = 20. Ini- 
tially, the density grows slowly all over the cluster volume, and this steadily enhances the 
radiative cooling, causing lower temperatures. Eventually, as the temperature approaches 
T~3 x 10 5 K, the cooling rate increases steeply, and thermal instability occurs. This low- 
ers the temperature rapidly to T = 10 4 K, particularly in the densest central regions. On 
the other hand, the outer regions, where the density is lower, remain hot and thus a large 
pressure gradient between the two regions leads to the formation of a strong shock wave 
propagating inwards (see Figure HI left column). 

The simulations show a very dynamic evolution in which regions of hot gas of different 
dimensions appear and grow close to the center, as more energy is deposited within the 
cluster volume. The hot gas expands super-sonically into the low pressure warm (10 4 K) 
surroundings. This decreases locally the density of the hot gas, preventing its becoming 
thermally unstable, while at the same time the parcels of warm gas are compressed from 
all sides into high density condensations, until they reach pressure equilibrium with the 
surrounding hot gas. The hot regions grow until they again occupy most of the cluster 
volume, see Figure HI middle column. 

While all this is happening, the inward propagating shock wave, overruns the central 
region, accelerating inwards most of the high density condensations, and these eventually 
leave the computational grid as they cross the inner boundary. As the high temperatures 
are reestablished, the overall pressure gradient vanishes and the shock wave decays. Once 
the hot gas again permeates the cluster volume, the density starts to grow again and the 
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whole process repeats itself, but in a weaker form because of the presence of a few dense 
condensations surviving from the previous cycle. This model exhibits 2- 3 subsequent weaker 



perio ds of oscillations similar to the ones observed in ID simulations (jTenorio-Tagle et al. 



20071 ; the low-energy model). 



Finally, the oscillations vanish and a quasi-equilibrium situation is established, in which 
high density condensations form at an approximately constant rate (see Figure HI right 
column). The amount of mass which goes into warm dense condensations is just enough to 
maintain the hot medium close to the thermally unstable regime. It is a quasi- stationary 
situation: the density tends to grow everywhere, surpassing some threshold value that favours 
the thermal instability in the densest central regions. This sudden loss of pressure prevents 
matter within the stagnation boundary from escaping the cluster as a wind. The thermally 
unstable gas is instead packed into high density condensations, most of which leave the 
computational grid through the inner boundary. The location of the stagnation surface, 
which separates the region where the thermal instability occurs from the outer stationary 
wind, also experiences some oscillations. This leads to a no n-spherical surface w ith an average 



radius close to that given by the analytic approximation (jWunsch et al.l 120071 ). 



d3 / T \ 1/2 



i? 3 

1 Lo 



sc \ -^sc 



1 - ■ (4) 



The R st radius at which vicinity the orientation of the velocity vectors abruptly changes from 
an outward to an inward motion has been indicated with a red line in the bottom panels of 
Figure HI 



3.3. Models 2-4, L sc /L crit = 2 

In these cases the stagnation radius Rst is smaller than in model 5, resulting in a smaller 
thermally unstable region. The lower mass deposition rate results in a smaller amount of high 
density gas being formed there (see Figure [5]). Otherwise, the evolution is similar to that of 
model 5, including the initial relaxation period of intense thermal instability, the appearance 
of re-pressurizing shocks which lead to high density condensations, and the exit of most of 
these through the inner boundary, to finally reach equilibrium between the formation of high 
density condensations and the mass deposition rate. 

We have computed this model on three different grids to check how the resolution affects 
the results (see the different rows of Figure [5]). Although the resultant fragments tend to be 
smaller and more structured on the higher resolution grids, the global characteristics such 
as the mass flux through the cluster border (see Tabled]) are in a reasonable agreement. 



- 15 - 



3.4. Model 6, L S c/L CTit = 200 

The hydrodynamical behavior of this case is very similar to that of model 5, the only 
differences are quantitative: In this case many more high density condensations form and 
occupy a larger fraction of the cluster volume (see Figure [6]). The quasi-stationary region, 
above the rapidly varying stagnation surface, becomes at times very narrow and, as it is 
repeatedly perturbed by high density condensations leaving the cluster, at times it is not 
even a contiguous region. Nevertheless, most of the thermally unstable gas driven into high 
density condensations stays within the cluster volume and after some time, as in model 5, 
goes through the inner boundary and disappears from the computational domain. 

3.5. Other Models 

The most important parameter of our cluster model is the ratio Lsc/^crit, as this defines 
the location of the stagnation surface and thus the relative sizes of the thermally unstable 
region and of the quasi-stationary outer wind region. We have performed many more sim- 
ulations (see Table 1), all presenting the same general features as in model 5. For example 
models 7 and 8, (L sc /L crit = 2.5, 77 = 0.1 and 0.3) have L sc /L crit values close to those of 
models 2-4. However, the smaller heating efficiency rj assumed in these cases, results in a 
smaller temperature and pressure of the hot medium, and this leads to a stationary wind 
that expands with a smaller velocity. The lower temperature of the hot gas leads to slower 
velocity re-pressurizing shocks. The lower ambient pressure of the hot gas also leads to larger 
final sizes of the high density condensations and after some time to their accumulation near 
the central zones of the computational grid. This sooner or later prevents the exit of matter 
through the inner boundary in case of models 9 and 10 (Lsc/^crit = 25, 77 = 0.1 and 0.3) 
and model 11 (L sc /L crit = 250, r\ = 0.1). We believe this is an artifact promoted by the fact 
that we do not allow for this matter to go into star formation. 

Models 12 and 13 are similar to models 5 and 9, respectively, but the gas is allowed to 
cool to 10 2 K instead of 10 4 K. A comparison between model 12 and 5 shows very similar 
results. This means that the reduced volume of clumps in model 12 (due to their lower 
temperature and hence pressure) does not affect the ratio of the clumps that are ejected 
from the cluster to the clumps that stay there and finally leave the computational domain 
through the inner boundary. Therefore, the outflow from the cluster (M out ) remains the 
same in both models. However, in the case of model 13 (Rsc — 3 pc, r] = 0.1), the reduction 
of the clump sizes prevents accumulating and filling the cluster with a dense warm material, 
as happens in model 9. As shown in Figure [7] and in Tabled] the formation of new clumps via 
thermal instability is compensated with the outflow from the computational zone making, 



model 13 more realistic. 
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Time = 0.025 Myr Time = 0.044 Myr Time = 0.400 Myr 
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Fig. 4.— Model 5 (L sc /L crit = 20, r] = 1) at t = 0.025 Myr (left column), t = 0.044 Myr 
(middle column) and t = 0.4 Myr (right column). The first two rows of panels show the 
logarithm of the wind temperature and density, respectively, across the whole computational 
domain. The third row shows the logarithm of pressure in the cluster central region (below 
Rsc — 10 pc), and the bottom row shows the wind velocity in the same region, as arrows 
together with its magnitude coded by the color scale. The red arc is the stagnation radius 
given by the semi-analytical solution. 
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Fig. 5.— Models 2-4 (L sc /L crit = 2, 17 = 1) at t = 0.25 Myr. The left and right columns 
show the logarithm of the wind temperature and density, respectively, in the cluster central 
region. The three rows of panels represent the different grid resolutions: the top row is 
150 x 56 (model 2), the middle row 300 x 112 (model 3) and the bottom row 600 x 224 
(model 4). 
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Fig. 6.— Model 6 (L sc /L CTit = 200, r] = 1) at t = 0.25 Myr. The top left and right panels 
show the logarithm of the wind temperature and density, respectively, across the whole 
computational domain. The bottom panels show only the cluster region (below i?sc — 10 pc), 
the left panel represents the logarithm of pressure, the right panel shows the wind velocity as 
arrows together with its magnitude coded by the color scale. The red arc is the stagnation 
radius given by the semi-analytical solution. 



-20- 




-1 -0.5 0.5 1 

r[pc] 



Time = 0.74 Myr 




-1 -0.5 0.5 1 

r[pc] 



Fig. 7. — Model 13 (Lsc/^crit = 25, r\ = 0.1) at t = 0.74 Myr. The meaning of the panels is 
the same as in Figure [61 
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Fig. 8. — The solid line represents the time evolution of the mass flux measured at the outer 
boundary of the computation domain in model 5. The mass deposition rate to the whole 
cluster M sc and its fraction deposited between R st and i?sc are represented by dash-dotted 
line and dashed line, respectively. The dotted line represents the average flux at the outer 
boundary in the period 0.2 - 0.8 Myr. Those average mass fluxes are shown as symbols in 
Figure M The rarefied wind produces the flux close to the value given by the dashed line. 
The maxima of M out are due to condensations passing through the outer boundary They are 
preceded by short periods in which M ont drops below the dashed line, as if the condensations 
"cast shadows" to the regions above them - the mass flux of the rarefied wind is slightly 
lower then and corresponds to the moments when the outer boundary is shadowed by an 
approaching condensation(s). 
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Fig. 9. — The fraction of the matter inserted within the stellar cluster which is able to leave 
it as a function of the cluster luminosity normalized by its threshold value (-ksc/Aait)- The 
solid line represents the fraction of mass which is inserted above R st and which leaves the 
cluster in a form of the hot wind (which may eventu ally cool down outside the cluster). The 
x symbols denote the ID simulations described in iTenorio-Tagle et al.l (120071 ). The other 
symbols represent the 2D simulations: triangles (models 2-6), the plus (model 7), the circle 
(model 8), the square (model 12), and the diamond (model 13). The values were obtained 
by averaging over time periods 0.2 — 0.8 Myr (models 2-6), and 1 — 2 Myr (models 7, 8, 
12 and 13). The dashed line is the fit to the outflows of 2D models. 
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4. Discussion 

Figure [8] shows the time evolution of the mass flux through the outer boundary of the 
computational domain for model 5. The outflow becomes quasi-stationary as an equilibrium 
is established between the formation of high density condensations and their ejection from 
the cluster and exit through the inner boundary. The average mass flux across the outer 
boundary of the computation domain (dotted line in Figure [8]) slightly exceeds the rate of 
mass deposition that occurs between the stagnation radius R st and the star cluster edge 
Rsc- This is because some dense condensations that formed in the thermally unstable inner 
region (r < R st ) cross the stagnation surface and eventually leave the cluster contributing 
to the total mass flux across the outer boundary of the computational grid. 

Figure [9] shows how the relative outflow from the cluster (measured as the mass flux at 
the outer boundary, averaged over long time periods) depends on the ratio Lsc/^crit- We 
plot all the models where the equilibrium between the clump formation and their removal 
through inner or outer boundaries is reached. Other models (Nos 9, 10 and 11), where the 
volume of the cluster is completely filled with the warm matter (see Sect. 3.5.) are omitted. 
The outflow from the cluster consists of two components: the hot wind, which originates in 
the outer region of the cluster above R st , and the warm dense clumps that are formed below 
Rst and that manage to pass to the outer region where they are ejected from the cluster. The 
figure shows three zones: the bottom zone is the fraction of the deposited mass that goes 
into the wind, the middle zone represents the mass in ejected high density condensations 
that stream away with the wind (< 20%Msc), and the upper zone is the mass which stays 
in the cluster (available for star formation). Note that the hot wind also cools down to 
temperatures T ~ 10 4 K at short distances from the cluster surface, so, finally there are 
also two phases in the outflow: the warm wind and the dense condensations (which should 
expand unless they cool even further). 

In order to advance the problem further, apart from an adequate hydrodynamic model 
that takes into consideration specific characteristics of SSCs (i.e. their high masses, small 
radii, large stellar densities and extreme output of mechanical energy), one would need to 
couple the hydrodynamics to the UV radiation field. We have assumed here that this is the 
case and that the UV radiation field generated by the massive stars evolving in the cluster 
keeps the temperature of the thermally unstable gas at T ~ 10 4 K through photoionization. 
This may be true for very young clusters, before the supernova era starts (tsN ~3 Myr). 
However, older clusters, with a reduced ionizing photon flux, would soon become unable 
to photoionize all the gas that became thermally unstable. In such cases, (see Figure [7j) 
the thermally unstable gas would continue to cool further, while being compressed into 
correspondingly smaller volumes. As a result, the increasingly high densities would trap the 
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ionization front in the outer skins of the condensations and their interior would remain neutral 
and at low temperature (~ 10 K). In this way, if a parcel of gas with an original temperature 
~ 10 7 K, that becomes thermally unstable and is able to cool down to 10 K, would experience 
a rapid evolution in which its volume, in order to preserve pressure equilibrium, would be 
reduced by six orders of magnitude while its density would become six orders of magnitude 
larger. 

Another factor not taken into account in the present set of simulations is gravity. The 
gravitational pull caused by the cluster is perhaps not significant for the high temperature 
gas, as this has a sound speed of several hundreds of km s" 1 , much larger that the escape speed 
from the cluster. However, it should promote a faster exit of low temperature condensations 
across the inner grid boundary. Indeed, if one considers a condensation that develops at the 
stagnation radius, its free-fall time to the cluster center will be r ff = vri?g/ 2 / (2GM st ) 1//2 where 
M st is the mass below R st . In pc and solar mass units, it is r ff « 16.5R^ /M^ Myr which 
leads to time-scales much shorter than the computational time. Thus gravity would lead to 
an increase in the speed of such condensations as they move to cross the inner boundary. The 
implication of this is a faster accumulation of the thermally unstable matter near the center 
of the cluster, where further generations of stars are to take place. Another important factor, 
also promoting a faster matter accumulation and further stellar generations arises from the 
self-gravity of the thermally unstable gas. 



5. Conclusions 

Here we have confirmed by means of 2D hydrodynamic simulations the existence of 
a bimodal solution for SSCs above the threshold line. We have shown that the evolution 
within the volume defined by the stagnation surface is very dynamic. The stagnation surface 
itself has a very dynamic morphology that continuously changes with time. Nevertheless, the 
average value of the stagnation radius remains close to the value predicted by ID simulations 
and semi-analytic solutions. This region suffers a very dynamic evolution in which parcels of 
gas continuously become thermally unstable and are rapidly driven into very small volumes to 
compensate their sudden loss of pressure. The number of such regions depends on the excess 
star cluster mechanical luminosity above the threshold value. The fraction of the cluster 
volume occupied by the warm medium depends on the balance between the formation of 
high density condensations, via thermal instability, and their removal via secondary star 
formation and/or their escape from the cluster. In our model, the secondary star formation 
is partly accounted for with the exit of high density gas across the inner boundary. However, 
a better treatment in the future would be to implement a more physical description of 



-25 - 



secondary star formation. 

We have also shown that the growth of high density condensations within a SSC volume, 
is strongly linked to the parameter rj, as this determines the location of the critical luminosity, 
L cr it, in the star cluster mechanical luminosity vs size plane, and thus determines how far 
above the critical luminosity a cluster is. It also influences the sound speed (and pressure) in 
the thermal instability region and as a result it fixes the final high density that condensations, 
confined by pressure, attain (higher t] =>■ higher pressure =>- denser condensations). 

We have shown that most of the condensations generated within the stagnation volume, 
are unable to leave the cluster volume and thus accumulate. Eventually this will result in 
further generations of star formation. The central implication of this result is that most of 
the metals processed by stars in massive and compact SSCs will not be ejected back into the 
host-galaxy ISM, an important issue to be taken into consideration by models of galactic 
chemical evolution and ACDM models of the universe. Careful analysis of observational data 
of the most massive and compact star clusters is now required to select those which may 
evolve in the bimodal, catastrophic cooling regime. For all of them we expect a mixture hot 
X-ray gas, a warm partially photo-ionized plasma, and an ensemble of accumulating cool 
dense condensations. Thus they should be detectable in the X-ray band, visible, radio and 
infrared regimes simultaneously. 

Note also, that the high stellar densities associated with SSCs resemble t hose of globular 



clust ers, in which star-to-star abundance inhomogeneities have been observed (IBekki fc Chiba 



20071 ). This can be understood if globular clusters have entered the bimodal regime during 
their early evolution, and this has led to multiple stellar generations forming with the matter 
reinserted into the cluster volume. Such clusters, having a reduced amount of matter being 
lost back into the ISM, would have remained more stable against disruption. 

Other open issues and some more astr ophysical con s equen ces of clusters undergoing this 



bimodal evolution have been discussed by lPalous et al.l (120081 ). 
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